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ABSTRACT 

The laminar-turbulent breakdown of a boundary-layer flow along a hollow 
cylinder at Mach 4.5 is investigated with large-eddy simulation. The subgrid 
scales are modeled dynamically, where the model coefficients are determined 
from the local resolved field. The behavior of the dynamic-model coefficients 
is investigated through both an a priori test with direct numerical simulation 
data for the same case and a complete large-eddy simulation. Both formula- 
tions proposed by Germano et al. and Lilly are used for the determination of 
unique coefficients for the dynamic model, and their results are compared and 
assessed. The behavior and the energy cascade of the subgrid-scale field struc- 
ture are investigated at various stages of the transition process. The investiga- 
tions are able to duplicate a high-speed transition phenomenon observed in 
experiments and explained only recently by the direct numerical simulations of 
Pruett and Zang, which is the appearance of "rope-like" waves. We also inves- 
tigate the nonlinear evolution and breakdown of the laminar boundary layer 
and the structure of the flow field during the transition process. 
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1. INTRODUCTION 

In high-speed boundary-layer flows, the physics, prediction, and control of the laminar- 
turbulent transition process are of increasing fundamental, as well as technological importance 
because of such national projects as the design of the High-Speed Civil Transport (HSCT) and 
the National Aero-Space Plane (NASP). In contrast to the vast progress that has been made in 
understanding the different mechanisms of laminar breakdown to turbulence in incompressible 
flows, progress has been much slower for compressible flows because of the inherent complex- 
ity of high-speed numerical and laboratory experiments. 

The road to transition in high-speed flows may involve several stages. First, an initial 
stage of linear instability of small disturbances, called primary instability, is well described by 
the compressible stability equations. (See Mack [1].) Multiple eigenvalues of high-frequency, 
amplified, and damped modes (Mack modes) appear in the solution of these equations. The 
first of the Mack modes is called the second mode, which is the most unstable of all the modes 
as a two-dimensional (2-D) disturbance. The linear stage is followed by a second stage, where 
apparent nonlinearity is revealed by the amplification of three dimensionality and the high 
growth rate of these disturbances. At this stage, the disturbances are still weak and the non- 
linear distortion to the mean flow can be negligible. Studies of secondary instability mechan- 
isms in high-speed flows (see El-Hady [2,3], Masad and Nayfeh [4], Ng and Erlebacher [5], 
Ng and Zang [6], and Pruett et al. [7]) have indicated that the growth rate of these three- 
dimensional (3-D) disturbances is more gradual; it persists for a long distance downstream, 
unlike the explosive and fast growth that was observed in incompressible flows. This may 
lead, in high-speed flows, to transition regions that are quite lengthy. These investigations 
have concluded that the subharmonic secondary instability is the most likely path to natural 
transition in high-speed flows. Finally, a strong nonlinear stage takes place with high and 
intense fluctuations that lead to the inevitable fully turbulent regime. 

Experimental studies at high speeds have been reviewed by Kendall [8]. These studies 
stressed the fact that the second mode is the dominant instability in the high-speed range and is 
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responsible for the production of persisting, well-ordered, "rope-like" waves seen in photo- 
graphs (e.g., Fischer and Weinstein [9]). In their direct numerical simulation (DNS) calcula- 
tions of a Mach 4.5 flow along a hollow cylinder at the beginning of the laminar breakdown 
stage, Pruett and Zang [10] were able to reproduce this rope-like structure. They numerically 
generated a schlieren photograph (like the laboratory experiment), which was derived from the 
wall-normal density gradient, and concluded that these flow structures are associated with the 
subharmonic secondary instability and arise from the 2-D projection of the staggered 3-D 
lambda vortices. 

The linear stability and secondary instability theories are not capable of describing the 
flow structure in the final stage that leads to laminar breakdown and afterwards. The recent 
approach of parabolized stability equations (PSE) [11,12] also has limitations and can describe 
the flow structure until just after the onset of transition. The flow structure during transition 
can only be described with DNS calculations. However, the application of compressible DNS 
to the laminar-turbulent transition process is an exceedingly expensive task because of the wide 
range of length scales present in the later stages of transition; hence, more resolution is needed 
in this stage than for laminar or fully turbulent flow. Typical temporal DNS computations 
require many hundreds of supercomputer hours; typical spatial DNS computations require 
thousands of hours. Because of its high computational requirement, DNS is usually used to 
simulate a forced transition to focus on the later stages of the transition process. 

With the increased interest in boundary-layer transition of high-speed flows, given the 
high cost of compressible DNS and the problems associated with high-speed laboratory experi- 
ments, other, less computationally intensive methods are needed for engineering purposes. A 
large number of ongoing efforts by different researchers are aimed at developing such less 
computationally intensive methods that will make a distinct technological contribution to 
expanding the state of the art in this field with large-eddy simulation (LES). 

The LES technique has been successfully applied to the study of turbulent flows (e.g., see 
the review article by Rogallo and Moin [13]) and has become an important method of 
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simulation. Ifi DNS, all flow scales are accurately resolved; in LES, only the large energy- 
carrying scales are accurately resolved. The influence of the discarded small or subgrid scales 
(SGS) must be modeled appropriately to simulate the energy cascade (i.e., the energy transfer 
between the large and small scales). A filtering process is usually used to separate the large- 
and small-scale motions, which results in the appearance of the SGS in the residual stress and 
heat flux terms. 

In LES, a unified SGS model is usually used for different types of flows because the 
small-scale motion is more universal for all types of flows. Specifically, in turbulent flows, the 
small-scale motions tend to be more isotropic and homogeneous; hence, they can be 
represented by simpler models, in comparison with the currently used Reynolds-averaged tur- 
bulence models. On the other hand, the large-scale motion differs sharply from one flow to 
another; it is directly affected by the boundary conditions. The main contribution of the model 
is that it allows the transfer of the correct amount of energy from the large to the subgrid 
scales, or vice versa, near the wall. 

A central issue in LES is the development of models for the small subgrid scales. Con- 
trary to turbulent flows, the SGS model used for transitional flows is required to capture 
several stages: the primary, the secondary, and the fully nonlinear stages. The model is also 
required to accurately predict the transition onset, the transition region, and the fully turbulent 
stage. A review of transition-region modeling was recently given by Singer [14]. The applica- 
tion of LES to transitional flows dictates numerous modifications and extensions for the con- 
ventional SGS models that are used for turbulent flows. Although the Smagorinsky eddy- 
viscosity model [15] remains the fundamental subgrid model in LES for incompressible flows, 
several extensions to this model have been made. Moin and Kim [16] used a damping func- 
tion to account for the near-wall effects. Piomelli and Zang [17] modified the Smagorinsky 
model by using a scale function to avoid the known excessive damping of the SGS motion. 
The study by Piomelli et al. [ 1 8] established the feasibility of using LES for transitional flows. 
They applied this technique to the incompressible flat-plate boundary layer with a modified 
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Smagorinsky eddy-viscosity model. The results compare quite well with the corresponding 
DNS calculations. 

The LES of high-speed flows has recently received increasing interest. For compressible 
isotropic turbulence, Speziale et al. [19] and Erlebacher et al. [20] have developed a compressi- 
ble analog of the linear combination SGS model of Bardina et al. [21] (called the SEZH 
model); Zang et al. [22] used this model and found good agreement with the DNS data for 3-D 
compressible turbulence. Comte et al. [23], Lesieur et al. [24], and Normand and Lesieur [25] 
have developed and used a structure-function SGS model that is based on the physical imple- 
mentation of the concept of eddy viscosity. They have applied the structure-function model to 
compressible isotropic turbulence, as well as to compressible transition. Krai and Zang [26] 
have tested both the SEZH model and the structure-function SGS model for supersonic tur- 
bulent flow over a flat plate and have concluded that the later model gives better agreement 
with experiments. 

In all previous studies, the ad hoc manner in which the model constants and the model 
modifications (to satisfy certain physical conditions) have been treated is not satisfactory, so 
that the question always remains: can we make LES as attractive as DNS by eliminating the 
need for ad hoc models? Recently, Germano [27] and Germano et al. [28] have answered this 
question by introducing a new SGS eddy-viscosity model (dynamic model) for incompressible 
flows that avoids the ad hoc treatment and utilizes the spectral information in the large-scale 
field to evaluate the small-scale stresses. They applied this model to both channel and 
boundary-layer flows. Moin et al. [29] have extended the dynamic model to compressible 
flows and scalar transport and have applied it to the LES of decaying isotropic turbulence. El- 
Hady et al. [30] have established the feasibility of using the dynamic SGS model in the LES of 
high-speed transitional flows. 

In this paper, we study the nonlinear evolution and breakdown of the laminar boundary 
layer in high-speed flows, as well as the structure of the flow field during the transition pro- 
cess. This study is achieved by using a large-eddy temporal simulation, where the large flow 
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structures are computed and the small scales are modeled. The flow field under study is the 
transitional axial flow field on an axisymmetric body (a hollow cylinder) at a Mach number of 
4.5 and at a Reynolds number of 10000, based on the boundary-layer displacement thickness. 
We also study some aspects of the dynamic modeling of the small-grid scales applied to this 
transitional boundary layer. We expect the interaction, or the energy exchange, between the 
large and the small scales to be quite different in compressible transitional boundary layers in 
comparison with previous applications to isotropic turbulence and incompressible transition. In 
the remainder of the paper, we introduce the governing equations, develop the SGS dynamic 
model, and present the numerical procedures for the LES. We then discuss the results of both 
the a priori test for the model and the LES for a transitional boundary layer at a Mach number 
of 4.5. Finally, we close with some concluding remarks. 

2. PROBLEM FORMULATION 


Governing Equations 

The axisymmetric flow field along a cylinder of radius R„ is described by the compressi- 
ble 3-D Navier-Stokes and energy equations written in the body-fitted coordinate system x,y, 
z, where x is the direction along the body, y is the normal to the body, and z is the azimuthal 
direction. Lengths, velocities, and time are made dimensionless with the boundary-layer dis- 
placement thickness 8* as a reference length, the edge velocity u‘, and 8*/« e *, respectively. The 
pressure is made dimensionless with p; u* 2 . The temperature, density, specific heats, viscosity, 
and thermal conductivity of air are made dimensionless with their corresponding edge values. 
The displacement thickness 8* is defined for the axisymmetric boundary layer by the quadratic 
equation 


8*[1 



* * * 



( 1 ) 


where r* is the radial coordinate defined as r*(y) = R* + y* . Equation (1) reduces to the 
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standard planar boundary-layer definition as R ‘ -> «. In terms of these dimensionless quanti 
ties, the governing equations in vectorial form read 

ie. 


dt 


+ V.(pV) = 0 


( 2 ) 


+ V.(pW) + A = -Wp + J-V.o + X 
af a 


+ V.Vp + yp V.V = — v.(* V0) + 

dt y IF M;TR R 


(3) 

(4) 


with the state equation for an ideal gas 

Y = p9 (5) 

In the above equations, a is the dimensionless viscous stress tensor, and <t> is the dimensionless 
dissipation function. They are defined as 


a . p[W + (W) r ] + X (V.V) I 

(6) 

<J> = o:W 

(7) 


Also, in this coordinate system, the vector operators are defined as 


w=^ + — ^ + - + - 

dx r dz dy r 


V.Vp . „#■ + + V& 


^ "dx r " dz dy 
and some extra terms A and % appear in the momentum equation; these are defined as 


( 8 ) 

(9) 


A = -^(0 , v , -w) (10) 

r 

( 11 ) 

Here, y is the ratio of specific heats (y = 1.4); M e is the edge Mach number; R = p e u e 8'/p e is 
the Reynolds number; r = c p \i e /k e is the Prandtl number (r = 0.72); p and X are the first and 
second coefficients of viscosity, respectively; k is the thermal conductivity; I is a unit tensor; 
and ( ) r denotes a transpose. The Stokes hypothesis is assumed, and the viscosity and thermal 
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conductivity are related to the temperature through Sutherland’s law. 

Instead of solving equations (2) through (5) directly for the full range of scales, we limit 
the resolvable scale size by explicitly applying the spatial filter G to each term in the Navier- 
Stokes and energy equations. This filtering process will separate the resolvable field from the 
SGS field, which is accomplished by decomposing each field quantity F in the flow domain D 
as 

F - F + F" (12) 


where the filtered quantity F is defined as 


F(x,t ) 


G(x-C> A)F(C,OrfC 


(13) 


and F" is the SGS part that accounts for those scales not resolved by the filter width A. We 
note here that the filtered quantity F in equation (12) or (13) is a spatially averaged, time- 
dependent resolved quantity, rather than a time-averaged or ensemble-averaged mean quantity 
as in the Reynolds-averaged Navier-Stokes equations. We keep A as a parameter in the 
integral (equation (13)) to clearly show that the structure of the large-scale field and, conse- 
quently, the structure of the SGS field will depend on the filter width. The smallest scale 
allowed by the filter and the SGS model, which is 0(A), should be greater than the smallest 
scale resolved by the grid size, which is 0(h). This relationship will ensure that the computa- 
tions are independent of the numerical algorithm used in the simulations. In these calculations, 
we choose the filter width A ; = 2 ft, , where ft, is the grid size in the i th direction. 

We use a sharp Fourier cutoff filter in these calculations and apply it in the wave space; it 
is conveniently defined as 

f i for ki < K ci 

Gj(k;) _ |^o otherwise (14) 


where Gi is the Fourier coefficient of the filter function in the / th direction and K ci is the cutoff 
wave number in the i th direction (related to the corresponding filter width A, by K ci = n/A ( ). 
For laminar-turbulent transition flow problems, a cutoff filter is more appropriate than other 
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commonly used filters (e.g., a Gaussian filter or a variable-width box filter). This filter ensures 
clear identification of the SGS with scales that have wave numbers larger than the cutoff wave 
number K ci , rather than a mixed identity with both small and large scales, which is the case for 
other filters (see Piomelli et al. [31]). 

To account for large density fluctuations in high-speed flows, the resolved velocity and 
temperature fields are written in terms of Favre-filtered quantities (Speziale et al. [19]), which 
are defined as 


F = pFIp 


(15) 


Now F is decomposed to 


F=F+F' (16) 

where F’ is the SGS part of F , based on Favre filtering; consequently, the Favre-filtered 
governing equations used for the large-eddy simulation are 

+ V.(pV) = 0 (17) 

+ v.(p^V) + A = -Vp + — V.3 - v.x + % (18) 

ot R 

+ V.Vp + to V.V = (Y - lX^Vp - V.Vp) + — l — V.(*W) + - ArV.q (19) 

dt * M?TR R M } 


with the state equation 


yM e 2 p = p9 


( 20 ) 


In deriving equations (17) through (20), we have utilized the commutative property of the 
filtering process with space and time such that 

VJ = V.F (21) 

dt dt 

Also, we have decomposed the stress tensor pVV into its resolved (pVV) and SGS components 
(x) and the heat flux into its resolved (pve) and SGS components (q). Also, in the filtering 
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process, the extra term A that results from the use of a body-fitted coordinate system will con- 
tribute to the SGS; it is now defined as 


A = -£*(0, x^,x H ) 


( 22 ) 


Equations (17) through (20) are further simplified by assuming that <3 = o(V), 
(y - 1XV.V/7 - V.Vp) = 0, k V9 = kVQ, O = 3>(V), and neglecting the contribution to the SGS by the 
viscous dissipation. 

The SGS stress tensor x and the SGS heat flux vector q appear in the right-hand side of 
the momentum and energy equations, respectively; they must be modeled to close these equa- 
tions. The SGS stress tensor x consists of the Leonard, cross, and Reynolds stresses, which are 
defined as 


Liu = P (M/ - ttk&i) 

Ckl ~ P ( u k % + &k u l ') 

Ru - p {u\u'i) 

The SGS heat flux vector q consists of similar components, which are defined as 

qL k = p (u*0 - u k Q) 
qC k = p (u\0 + Ute 7 ) 
qR k = p (u 9 #) 


(23) 

(24) 

(25) 

(26) 

(27) 

(28) 


Supplemented with appropriate initial and boundary conditions, equations (17) through (20) are 
used to yield the resolved flow field for all later times. 


Modeling the Small-Grid Scales 

To model the SGS stress tensor x and the SGS heat flux vector q, we use the approach 
given by Speziale et al. [19] and Erlebacher et al. [20], which is based on the eddy- viscosity 
Smagorinsky model [15] and an analogous eddy diffusivity model; these are modeled as 
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Xu = P (Ml -3,5,)- 2C r pA 2 l S\S kl '+ |c, pA 2 l S 1 2 6„ 

(29) 

q k =p(« t 0-5,e)-C e pA 2 ISl^- 

OXk 

(30) 

where 5„ is the Kronecker delta. The characteristic filter width A and the Favre-filtered rate- 
of-strain tensor S u are defined as 

A = (A, A y A z ) m 

(31) 

~ du 

“ " 3* 

(32a) 

= l*i 3« 

* " 2 r dz dx } 

(32b) 

2 { dy dx } 

(32c) 

= R a d\V , V 

r dz r 

(3 2d) 

~ 1 3v 3w n> 

^ “ 2 ( r dz + dy t 

(32e) 

5 = -^- 

w a>’ 

(320 

and 


Sid ' = S*/ “ 

(33) 

is 1 = (2 S ki Su ) 1/2 

(34) 

Note that the first term in equations (29) and (30) represents the scale similarity part of the 
model; it is the sum of the Leonard and cross portions of the SGS stress or heat flux fields. 
The remaining part of the equation represents the model for the Reynolds portion of the stress 
or heat flux fields. Piomelli et al. [31] indicated the need for the model to be consistent with 
the type of filter used in the analysis. With the cutoff filter, the scale similarity portion of the 
model is neglected, and the SGS field is modeled by the eddy-viscosity and its analogous 
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eddy-diffusivity portions. The coefficients of the model C r , C,, and C e are to be determined. 

The Dynamic Model 

In accordance with Germano et al. [28] and Moin et al. [29], the filtered equations of 
motion (equations (17) through (20)) are filtered once more. Again, a sharp Fourier cutoff 
filter is applied in the wave space (see equation (14)), where A; represents the new filter width 
in the i th direction and the ratio 5 = A,/A ; > 1. A ratio of 5 = 2 is adopted in the present calcu- 
lations. The filter is applied in the * -direction along the body and in the z -azimuthal direction; 
no filter is applied in the direction normal to the wall. The first and the second filters, usually 
called grid and test filters, respectively, will each produce a resolved flow field. The difference 
between the two resolved fields is the contribution by the small scales of length between the 
first and the second filter widths. In deriving a dynamic model, several assumptions are made, 
which we will identify. First, the same type of filter (a sharp Fourier cutoff filter) is used in 
both the grid and test filtering. Second, the same functional form, based on the Smagorinsky 
eddy-viscosity closure, is used to model the SGS stress terms from both filters. Third, the 
same functional form based on an analogous eddy-diffusivity closure is used to model the SGS 
heat flux terms from both filters. Fourth, a plane-averaging process is used in the derivation of 
the model coefficients, where the averaging is taken over a plane that is parallel to the wall. 
Finally, the model coefficients are independent of the filtering process. Each of these assump- 
tions has its merits and can open the way for future adjustments of the dynamic model (e.g., 
Ronchi et al. [32] and Ghosal et al. [33]). Here, we briefly describe the derivation of the 
model and leave the details to other references. (See references [28] and [29]). 

In accordance with Germano et al. [28] and Moin et al. [29], we apply a test filter to the 
equations of motion (17) through (20). We again use a sharp Fourier cutoff filter in the x- 
direction along the body and in the s -azimuthal direction (with no filtering in the y -direction 
that is normal to the wall). This new filter has the characteristic width 


A = (8a* A, 8Aj ,/? 


(35) 
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The SGS stresses and heat flux that appear in equations (18) and (19) can be defined as 


% = p (k*k, - u k u,) 


Qk = P - «*§) 


( 36 ) 


(37) 


Analogous to the above equations, the SGS stresses and heat flux that appear in the test-filtered 
equations can be defined as 

(38) 

(39) 


. * * 
T u = P («t«< “ «*“/) 


/V 


Qk = P («*0 - K*0) 


By assuming the same functional form and model coefficients and neglecting the similarity 
portion, the "test" SGS stresses and heat flux are modeled the same way as in equations (29) 
and (30), with the scale similarity portion neglected as 


T u = -2C r pA 2 l §\§ u ’+ jC, pA 2 l §\ 2 b u 

Qk =-C B pA 2 \§\^- 


(40) 

(41) 


The grid-filtered SGS field and the test-filtered SGS field are related by the Germano identity 
[27], such that 


/ — A A 

D k , £ T u - t u = p u k u, - piikpui/p 
E k = Q k - q k = pu k Q - pu k pQ/p 


(42) 

(43) 


where D u and E k are computable from the resolved large-eddy field. They are the resolved 
components of the stress tensor with scales of motion between the test and the grid scales. 
Lilly [34] called it the test window. 

By using the traces of equations (29), (40), and (42), the following expression can be 
derived to solve for C, 

= 2 C, (A 2 p I § 1 2 - AW? 2 ) = C,B (44) 


By substituting equations (29), (40), and (44) into equation (42) and equations (30) and (41) 
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into equation (43), the following expressions can be derived to solve for C r and C 0 : 

Du — 8« = — 2C r (A 2 ^ I S I Su — A 2 pl S I Sy ) = C r Afu (45) 

E k = -C e (A 2 $l S I - A 2 j5l 5 I -|^-) = C e W* (46) 

axjt ox k 

Each of the expressions (45) and (46) represents a set of six independent equations in one unk- 
nown. In order to uniquely determine the model coefficients C, and C 0 , we use the contrac- 
tions proposed by Germano et al. [28] and Moin et al. [29] in which equation (45) is con- 
tracted with Su and equation (46) is contracted with 3§/3x*. We refer to this contraction as c i 
in our calculations. We also use another contraction obtained from the least-squares analysis 
of Lilly [34] in which equation (45) is contracted with M u and equation (46) is contracted with 
N k . We refer to this contraction as c 2 in our calculations. The model coefficients in expres- 
sions (44), (45), and (46) are functions of space and time. They can locally become negative, 
which allows for backscatter (i.e., the transfer of energy from the SGS to the large scale). In 
their DNS calculations of transitional and turbulent channel flows, Piomelli et al. [35] have 
shown that about half of the grid points experience a backscatter when a Fourier cutoff filter is 
used. To alleviate an observed numerical instability caused by ill-conditioned local quantities, 
we perform a spatial averaging technique in a plane parallel to the wall (indicated by < >). 
This plane averaging enhances the numerical stability of the simulations at the expense of los- 
ing some of the conceptual advantages of the dynamic modeling formulation. The plane 
averaging leads to the following expressions for the model coefficients: 


c ( . <D “ > 

(47) 

< B > 


< (Du - D w 8«/3) %ki > 

(48) 

< M u Xu > 


„ < D k A k > 

C a — m T , 

(49) 

<N k A k > 



where Xu and A k assume the value given by the type of contraction. With the calculation of 
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the model coefficients, equations (29) and (30) are used to calculate the SGS stresses and the 
heat flux to close the governing equations (17) through (20). 

3. TESTING THE DYNAMIC MODEL 

In this section, we investigate the behavior of the dynamic model coefficients using the 
DNS data of Pruett and Zang [10] for a transitional boundary-layer flow along a hollow 
cylinder at a Mach number of 4.5. This so called a priori test serves as a qualitative evaluation 
of the model coefficients, the modeled SGS stresses, and the heat flux at various stages of the 
laminar-turbulent transition process. It also serves as an assessment of the two types of con- 
tractions that are used in calculating a unique value for the coefficients of the dynamic model. 
But, one should note that the output of an a priori test may not include the anticipated interac- 
tion between the modeled scales and the large scales, which can only be achieved through a 
complete LES. 

The initial resolution of the DNS data at time t = 0 was 12 x 6 x 96 in the stream wise, 
azimuthal, and wall normal directions respectively (symmetry was enforced about the plane z = 
0 to reduce the number of grid points by half in the azimuthal direction). As time evolved, the 
streamwise and azimuthal grid refinements were made as necessary to maintain at least a 
difference of eight orders of magnitude between the most and the least energetic Fourier har- 
monics. A well-resolved flow field was reached until the period 55 (time is given here in terms 
of the period of the primary wave) with a resolution of 96 x 48 x 144, beyond which the simu- 
lation became increasingly ill-resolved. 

We use the DNS data to examine the model coefficients on a coarse grid that is overlayed 
on top of the fine DNS grid. The effect of the filter width is investigated first by applying the 
grid filters 48 x 24 x 96, 32 x 16 x 96, and 24 x 12 x 96 to the DNS data (which correspond to 
Fourier cutoff wave numbers of 24, 16, and 12, respectively, in the streamwise and azimuthal 
directions). To evaluate the dynamic model coefficients, test filters are then applied to the 
DNS data and the ratio 8 = 2 is kept fixed. All calculations in this paper enforce a symmetry 
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about the plane z = 0 to reduce by half the number of grid points in the azimuthal direction. 
The numerical grid is always given in this order: streamwise, symmetric azimuthal, and wall- 
normal directions. 

We find that the model coefficients C r , C ; , and C e have a zero value at the wall and van- 
ish in the free stream; they oscillate at locations in the wall region. The coefficients are very 
small in the early linear region, and their values increase sharply as the laminar breakdown 
begins. Figure 1 shows the effect of the filter width on the plane-averaged values of the model 
coefficients at / = 54 inside the transition region. Their values increase as the filter width 
increases (coarser grid). A very coarse grid results in a sharp rise in the value of the 
coefficient with high oscillations in the wall region. This result will have a direct effect on the 
amount of energy contained in the SGS as the filter width changes. For example. Fig. 2 shows 
this effect on the plane-averaged values of the component x xl of the SGS Reynolds stress ten- 
sor, the SGS kinetic energy K, and the component q z of the SGS heat flux vector. Contraction 
ci is used in the calculations of Figs. 1 and 2. Figure 2 shows that the SGS stress, heat flux, 
and energy are very small for the grid 48 x 24 x 96 (as compared to the corresponding 
resolved quantities in Fig. 17 below), which indicates that with the Fourier cutoff filter the grid 
48 x 24 x 96 might be capable of resolving the flow field without modeling. 

The ratio 8 is the only input parameter needed for the dynamic model. For a fixed-grid 
filter, the effect of changing 8 on the model coefficients is expected to be the same as chang- 
ing the grid filter in Fig 1. As 8 increases, the test-filter width increases which results in high 
oscillations in the model coefficients and, consequently, in the SGS field. Figures 3 and 4 
show that effect on the model coefficients and the SGS field, respectively, at t = 54 using con- 
traction c 2 . The ratio 8=3 gives high oscillations in the wall-normal distribution of C r , C 9 , and 
the SGS Reynolds stress and heat flux. A ratio of 8= 1.3 results in negative values of C, and 
the SGS turbulent kinetic energy. The ratio 8 = 2 is kept fixed throughout all the calculations. 

An example of the effect of the proposed contractions c x and c 2 is shown in Figs. 5 and 6. 
Figure 5 illustrates the effect on the coefficients C r and C e . Figure 6 shows the effect on the 
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component x xz of the SGS stress and on the component q t of the SGS heat flux, which is given 
for a grid filter 32 x 16 x 96 at t = 54. The effect is a typical example for other grid filters and 
for different times. The contraction c 2 obtained from the least-squares analysis produces model 
coefficients, as well as SGS stress and heat flux components, that are more than 50 percent less 
than their counterparts with the contraction c x . Note that the model coefficient C, and, there- 
fore, the SGS kinetic energy are not affected by the type of contraction in an a priori test. We 
will show in section 5 that this may not be the case for a complete numerical simulation. 
Although the numerator and the denominator of the expressions in equations (47) through (49) 
are dealiased, the high oscillations persist, as shown in the figures for the model coefficients 
and the SGS field. These oscillations are probably unphysical and caused by ill-conditioned 
behavior. 

4. NUMERICAL PROCEDURES FOR LARGE-EDDY SIMULATIONS 

The solution of equations (17) through (20) uses basically the same algorithm used by 
Erlebacher and Hussaini [36] to simulate compressible flow over a flat plate, and by Pruett and 
Zang [10] to simulate compressible flow over axisymmetric bodies. Here, we briefly describe 
the main features of the algorithm. 

Because the parallel mean flow does not satisfy the governing equations (17) through 
(20), forcing terms are added to the right-hand side of these equations to suppress the temporal 
evolution of the unperturbed mean flow, so that the laminar state becomes stationary. The ini- 
tial value problem defined by equations (17) through (20) is explicitly integrated in time by 
using a third-order low-storage Runge-Kutta method. The time step is automatically controlled 
by the minimum of the maximum time steps allowed by the linearized advection and diffusion 
problems. The streamwise and the azimuthal directions arc periodic so that the flow variables 
can be represented by Fourier series in these directions. The periodicity lengths in the stream- 
wise and azimuthal directions (2 L x and L z ) are related to the corresponding wave numbers a 
and p by L x = 2n/a and L, = 2n/p. The spatial derivatives in these directions are evaluated at 
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each time step by spectral collocation methods. In the wall-normal direction, a sixth-order 
compact difference scheme is implemented. To increase the resolution requirements in the 
direction normal to the wall without drastically decreasing the time step in an explicit scheme, 
Pruett and Zang [10] have implemented a sixth-order compact difference scheme in the wall- 
normal direction, instead of the Chebychev collocation method used by Erlebacher and Hus- 
saini [36]. 

The grid is uniform in the x- and z -directions; a mapping is used in the y -direction, which 
clusters points near the wall and the critical layer and stretches the grid toward the far-field 
boundary. Symmetry is enforced about the plane z = 0 to save computational time. Spectral 
methods are known to introduce both truncation and aliasing errors. To minimize these errors, 
we adopted the same method as in Ref. [10]; we check the decay of the spectrum of the 
Fourier coefficients and refine the grid accordingly as needed. 

Equations (47) through (49) are used to evaluate the model coefficients at each time step 
from the flow field at the previous time step. The SGS stresses and heat flux are evaluated 
with equations (29) and (30) at each time step and added to the right-hand side before the 
governing equations are advanced in time. 

At each time step, the energy content E aP in any Fourier mode (a,p), normalized with the 
mean-flow energy E 0 is calculated from the Fourier coefficients of the velocity v a P as 

£a,u ^ = "iclf | POO I Vp(y4)l 2 ^jp~ d y 

<4=2- 8*o (50) 

Simulations are carried out for the flow along a cylinder with a Mach number of 4.5, a stream 
temperature of 110°R at the edge of the boundary layer, and at a streamwise station that 
corresponds to R = 10000. 

The initial conditions for the temporal numerical simulation match those of the DNS of 
reference [10]. They consist of the mean flow over a hollow cylinder q 0 (y), an axisymmetric 
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second-mode primary disturbance q^xys), and a secondary subharmonic disturbance q 2 {x,y,z,t). 
The mean flow is a spectral solution of the axisymmetric compressible boundary layer [37]. 
The solution of the primary disturbance [5] assumes the form 

q { (x,y,t) = 6] [^,(y> e' (ax ~ m) + cc] (51) 

and the solution of the secondary disturbance [7] assumes the form 

q 2 (x,y,z,t) = e 2 [* ia( ^ ,a)/2 e m ° Z £ 9 2 *(y) + cc] (52) 

n ~-N 

where a = 0.2523 and w = 0.2285 + i 0.319 x l(T 3 are, respectively, the real wave number and the 
complex frequency of a primary disturbance of amplitude ej = 0.085. The secondary disturbance 
has four wave components (N = 4), an azimuthal wave number fl = 0.2673, and an amplitude 
e 2 = 0.017. The eigenfunctions of both the primary and secondary disturbances are normalized 
such that their amplitudes measure the maximum root- mean-square value of the corresponding 
temperature disturbance. 

Dirichlet boundary conditions are applied to the total flow (except density) at the wall and 
at the far field. The density is calculated by projecting forward in time with the Runge-Kutta 
integrator, and the state equation is used to calculate the pressure. 

5. RESULTS OF LARGE-EDDY SIMULATIONS 

Features of the Dynamic Model 

A large -eddy simulation is conducted with the dynamic model and the initial grid 24 x 12 
x 96 by using the contraction c 2 and 8=2. The grid is fixed at 24 x 12 x 96 throughout the 
transition region. Figure 7 shows the time evolution of the dynamic model coefficients C r , C,, 
and C e . Although the coefficients are functions of y and /, they are averaged in the wall- 
normal direction for the purpose of this graph. The figure illustrates one of the favorable 
characteristics of the dynamic model; it turns itself off and on without the need for any added 
ad hoc functions. The model is turned off in the linear and weakly nonlinear regions. These 
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regions are approximated from the evolution of the primary and secondary disturbances; the 
linear region is between the period 0 and the period 15, and the weakly nonlinear region fol- 
lows until the period 35. The model turns itself on in the strong nonlinear region (from 
periods 35 to 45). The value of the model coefficients gradually increases as the laminar break- 
down commences, continues to increase through the transition region, and then decreases again 
in the fully turbulent region. 

A second LES is conducted in which a grid of 24 x 12 x 96 is used at the initial stages, 
and is refined as necessary to maintain at least a difference of four orders of magnitude 
between the most and the least energetic Fourier harmonics. The final grid reached is 48 x 24 x 
96. The time evolution of the dynamic model coefficients C r , C,, and C e (averaged in y) is 
included in Fig. 7 for comparison with the first LES of grid 24 x 12 x 96. 

For the LES with final grid 48 x 24 x 96, we examine the evolution of the wall-normal 
distribution of the model coefficients with time. Figure 8 shows the y -variation of plane- 
averaged values of C f , C,, and C e at t = 45, t = 54, and t = 62, which correspond to nearly the 
beginning, the middle, and the end of the transition region. First, the figure shows that the 
model gives the proper asymptotic behavior near the wall and vanishes in the free stream 
(without the need for ad hoc damping functions). As y approaches the free stream, we notice 
in the calculations that either the numerator of che formula for the model coefficient goes to 
zero or both the numerator and the denominator go to zero. In the later case, we assign a zero 
value for the model coefficient. Second, the model is first turned off and then turned on 
automatically in the transition region (without the need for an ad hoc intermittency function). 
Third, the values of the model coefficients are almost negligible in the linear and early non- 
linear stages and increase sharply in the wall region as the breakdown progresses. As the tur- 
bulent stage approaches, the values of the model coefficients drop and are not confined to the 
wall region, but spread away from the wall (not shown). At any plane parallel to the wall, the 
model coefficient C,{x,z,t) is always positive by the definition of the SGS kinetic energy, 
although negative values of C, sometimes occur during the calculations (see Fig. 8) in the 
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linear and the early nonlinear stages when the SGS kinetic energy is very small. This result 
leads to speculation on the validity of using the dynamic model to accurately represent the 
physics of flow fields that are dominated by the large-scale structure. On the other hand, model 
coefficients C r (x,z,t ) and C e (jt,z,/) can take a positive or negative value, but their plane- 
averaged values are rarely negative. 

The plane-averaged value of the turbulent Prandtl number Pr, is defined by the ratio C r /C e . 
This parameter directly influences the modeling of the SGS heat flux in Eq. (30). In the y- 
region, where both C r and C e have a nonzero value. Fig. 9 shows the evolution of the wall- 
normal distribution of the turbulent Prandtl number with time. The figure indicates that Pr, 
reaches an average value of 0.2 to 0.5 in the transition region. 

Contraction c x is used in the simulations that produced Fig. 8. When contraction c 2 is 
used, similar observations are noted in addition to the reduction in the value of the model 
coefficients predicted by the a priori test. Figures 10 and 11 show the effect of the contraction 
type on the model coefficients at t = 54 and t = 62, respectively. Although the type of contrac- 
tion does not affect the value of C, in an a priori test, a large-eddy simulation does affect their 
values, most probably because of the interaction between the modeled and large scales. 
Although the values for C r and C e predicted by the contraction c x are higher than those 
predicted by the contraction c 2 , the value of C ; calculated with contraction c x is slightly lower 
than that calculated with contraction c 2 as shown in Figs. 10 and 11. 

Global Features of the Flow Field 

For the LES with a fixed grid at 24 x 12 x 96 throughout the transition region, we 
checked the one-dimensional energy spectra against the corresponding spectra of a coarse-grid 
DNS (no SGS model) and a fine-grid DNS. The coarse-grid DNS seems capable of resolving 
the flow field only to a stage near the onset of transition. Its one -dimensional energy spectra in 
both the streamwise and the spanwise directions gives good agreement with the fine-grid DNS 
until t = 50. Beyond t = 50 and inside the transition region, the disagreement is slightly 
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enhanced by the dynamic SGS model because of the large increase in the value of the model 
coefficients. However, the dynamic model with the 24 x 12 x 96 grid seems unable to provide 
the right amount of energy transfer between the large and subgrid scales, as shown in Fig. 12 
at t = 55. Also, the energy content of some of the principal modes (Fourier harmonics) is 
checked against the corresponding energy content of a fine-grid DNS. This check demonstrates 
that the dynamic model with the grid 24 x 12 x 96 is capable of predicting the correct energy 
content of various Fourier modes only to a stage near the onset of transition. 

For the second LES calculations, where we reached a final grid of 48 x 24 x 96, Fig. 13 
shows the one-dimensional spectra of the kinetic energy in both the streamwise and spanwise 
directions at t = 55 (where the grid is 36 x 24 x 96). The figure exhibits good agreement with 
the fine-grid DNS. 

Figure 14 shows the time evolution of the principal Fourier harmonics in terms of its 
energy content. Of these, three harmonics ((1,0), (1/2,1), and (3/2,1)) are present in the initial 
conditions (together with their complex conjugates); however, the others ((0,2), (1,2), and 
(2,2)) emerged from the self-interaction of different components of the secondary disturbance. 
The different stages of evolution (linear, weakly nonlinear, strongly nonlinear, and breakdown) 
are recognized in accordance with the development of the primary mode (1,0) and the secon- 
dary subharmonics (1/2,1) and (3/2,1). The evolution of the various modes agrees well with 
the fine-grid DNS [10]. Of interest is the LES prediction of the domination of the mode (0,2) 
near the period 30 as indicated by the DNS [10]. 

In Figs. 15 and 16, we show a close comparison inside the transition region between the 
two LES presented in this paper (with grids 24 x 12 x 96 and 48 x 24 x 96) and the fine-grid 
DNS 96 x 48 x 144 of Pruett and Zang [10] for the energy content of the principal modes 
(1,0), (1/2,1), (0,1), and (0,2). The dynamic model with the grid 48 x 24 x 96 shows good 
agreement with the DNS and is capable of predicting the correct energy content of various 
Fourier modes all the way to and through the transition region; as mentioned before, in the 
case of the dynamic model with the grid 24 x 12 x 96, the comparison stops near the onset of 
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transition. 

The evolution of the model coefficients will directly influence the development of the 
SGS stresses and the SGS heat flux. However, the local values of the SGS field are not essen- 
tial to the evolution of the large-scale field. The divergence of the SGS field is the quantity 
that enters the momentum and energy equations. Figure 17 shows the development of the 
resolved field with time. The figure gives the y -variation of the plane-averaged values of the 
Reynolds stress component x xz , the kinetic energy K, and the heat flux component q z at pro- 
gressive times. The figure also shows a comparison with the fine-grid DNS results at t = 45 
and t = 54, where the DNS results are filtered on the same grid that is used for the LES. 
Fine-grid DNS results are not available at t - 62. 

Figure 18 shows the effect of the contraction type on the resolved field and a comparison 
with the fine-grid DNS results at t - 54. The fine-grid DNS results are filtered on the same 
grid that is used for the LES. Although both types of contraction resulted in a resolved field 
that is slightly overpredicted at this time ( t = 54), contraction c 2 is better able to predict the 
peak values than contraction c,. 

The time evolution of the shape factor H and the plane-averaged skin friction C f are 
shown in Fig. 19. The results from both types of contraction are compared with a fine-grid 
DNS and a coarse-grid DNS. The coarse-grid DNS is the same grid used for LES, but without 
the SGS modeling. A remarkable agreement exists between LES calculations and the fine 
DNS results, which indicates that the dynamic model ensures the transfer of the correct amount 
of energy between the large and the subgrid scales. Without the SGS dynamic modeling, the 
coarse-grid DNS results fail to predict the evolution of the skin friction during the transition 
region. Notice that both types of contraction used in the LES calculations give the same 
results until the end of the transition region, where some differences appear. 
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The Structure of the Flow Field 

The subharmonic secondary instability is the most likely path to natural transition in 
high-speed flows [2-7]. This instability leads to a flow-field structure that is dominated by 
periodic lambda-shaped vortices that are shifted by one half of a wavelength in the spanwise 
direction and staggered in the streamwise direction. Figure 20 shows a plan view, just below 
the critical layer, of the total vorticity predicted by the LES with the contraction c\ and by the 
fine-grid DNS at t = 45. The figure indicates an excellent prediction of the flow structure by 
LES at this time. The contraction c 2 gives identical results (not shown). 

The rope-like wave structure that was observed in high-speed experiments and was 
explained recently by DNS is associated with the subharmonic secondary instability [10] that is 
manifested in these computations at the time t = 45. This structure is the result of a 2-D pro- 
jection of the 3-D lambda vortices. The LES is able to capture the rope-like wave structure, 
even with the coarse grid (24 x 12 x 96). Figure 21 compares the rope-like structure that is 
derived from the wall-normal density gradient at t = 45 with the corresponding structure from 
the fine-grid DNS of Pruett and Zang [10] with the 64 x 36 x 128 grid. In Fig. 21 and in the 
following comparisons with DNS, the DNS grid is filtered to the corresponding LES grid. The 
significance of Fig. 21 is that although a coarse-grid DNS (24 x 12 x 96, which is the same as 
the LES grid) might be able to capture this rope-like structure, LES, with the dynamic model, 
causes no damping to the flow structure during the early stages of transition as does the stan- 
dard Smagorinsky model [15]. 

Following the periodic lambda-shaped vortices, detached high-shear layers of the stream- 
wise velocity are subsequently developed, which ride on top of the lambda vortices. With the 
intensity of the high-shear layers increasing, they stretch and eventually roll up. Sharp, low- 
velocity pulses (spikes) are then observed in the peak positions. At several stages in the transi- 
tion region ( t = 50, t = 53, t = 54, and t = 55), we compare the flow structure, in terms of 
the spanwise vorticity that results from LES calculations, with the DNS results [10], Figures 
22-25 show this comparison for both contractions used in the LES calculations. These figures 
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illustrate several points. First, excellent agreement exists between the LES and DNS results; the 
LES, with the dynamic model, is able to capture the bulk of the flow structure. Second, both 
contractions give identical results at the early stages of the transition region (at t = 50, the 
maximum spanwise vorticity is 4.9 for both contractions, compared with 5.1 for DNS) Third, 
slight differences in the details of the flow structure begin to occur in the late stages of transi- 
tion. Fourth, at these late stages of transition, the LES flow-field structure with contraction cj 
compares more favorably with the DNS results. 

As the transition process continues, small structures are formed due to the breakdown of 
the high-shear layer. This results in the periodical appearance of local regions of turbulence in 
the streamwise and spanwise directions. These regions spread as they travel downstream until 
a fully turbulent flow is attained. The later developments of the flow field are never under- 
stood in high-speed flows because of rare DNS calculations that cover the high computer 
demand of this stage. In Fig. 26, we present the results of the LES with contraction c t ; these 
results describe the spanwise vorticity structure of the flow field during the entire transition 
region and for the beginning of the turbulent region. The figure shows different stages: the for- 
mation of the lambda vortices, the development of the high-shear layers, the appearance of 
spikes, the decay of the large-scale structure, and the formation of new, small structures. 

6. CONCLUDING REMARKS 

The dynamic eddy-viscosity subgrid-scale model of Germano [27] and Germano et al. 
[28] has been applied to a high-speed transitional boundary layer at a Mach number of 4.5. In 
this method, the model coefficients C r , C, , and C 0 are determined dynamically as the calcula- 
tions proceed. The values of the model coefficients are adjusted automatically at different time 
steps with the local information contained in the resolved flow field. 

The application of the dynamic model to a high-speed transitional boundary layer is suc- 
cessful. The model gives the proper asymptotic behavior of the modeled quantities near the 
wall and in the free stream. The model has no dissipative character like the standard 
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Smagorinsky model; on the contrary, it turns dynamically off in the linear and early nonlinear 
regions and on both as the flow approaches the breakdown and in the transition region. The 
LES with the dynamic model is able to capture the "rope-like" wave structure, even with the 
coarse grid 24 x 12 x 96. 

The LES with the dynamic model is able to capture the bulk of the flow-field structure. 
While both contractions give identical results at the early stages of the transition region, slight 
differences in the details of the flow structure occur in the late stages of transition. For these 
stages, the LES flow-field structure with contraction c x compares more favorably with the DNS 
results. 

Large-eddy simulation of transitional flow along a cylinder at Mach 4.5 is achieved with 
one sixth of the grid resolution that was used for DNS. The structure of the flow field during 
the entire transition region and the beginning of the turbulent region is demonstrated with LES. 
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Figure 1. Effect of filter width on plane-averaged coefficients of dynamic model at t 
(inside transition region) with contraction c ^ (a) C r , (b) C, , and (c) C e . 
























Figure 5. Effect of different contractions on plane-averaged coefficients of dynamic model for 


grid filter 32 x 16 x 96 at t = 54: (a) C r , (b) C e . 





Figure 6. Effect of different contractions on plane-averaged (a) component x a of the SGS stress 
tensor and (b) component q, of heat flux vector for grid filter 32 x 16 x 96 at t = 54. 
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Figure 7. Evolution in time of plane-averaged coefficients of dynamic model with contraction 
c 2 : (a) C r , (b) c,-, and (c) C 0 . Coefficients are averaged in wall-normal direction for this graph. 
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Figure 12. One-dimensional spectra of flow kinetic energy at t = 55 for 24 x 12 x 96 LES with 
contraction c x : (a) streamwise spectra and (b) spanwise spectra. 






Figure 13. One-dimensional spectra of flow kinetic energy at t = 55 for 48 x 24 x 96 LES with 
contraction (a) streamwise spectra and (b) spanwise spectra. 
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Figure 14. Evolution in time of principal Fourier harmonics in terms of energy content for LES 
with final grid 48 x 24 x 96 and contraction c j. 
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Figure 15. Comparison of time evolution of energy content of Fourier harmonics (1,0), (1/2,1) 
inside transition region between LES (24 x 12 x 96), LES (48 x 24 x 96), with contraction C[ 
and fine-grid DNS (96 x 48 x 144) [10]. 
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Figure 16. Comparison of time evolution of energy content of Fourier harmonics (0,1), and 
(0,2) inside the transition region between LES (24 x 12 x 96), LES (48 x 24 x 96), with con- 
traction ci and fine-grid DNS (96 x 48 x 144) [10]. 
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Figure 17. Evolution with time of large-scale field (resolved) with contraction c, in terms of 
plane-averaged (a) Reynolds stress component x a , (b) kinetic energy K, and (c) heat flux com- 


ponent q, . 
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Figure 18. Comparison of large scale results (resolved) at / = 54 with contractions c\ and c 2 
and fine-grid DNS results: (a) Reynolds stress component (b) kinetic energy K, and (c) heat 
flux component q , . 
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Figure 19. Comparison of time development of (a) the shape factor H and (b) the skin friction 
Ct with contractions cj and c 2 with both fine- and coarse-grid DNS results. 
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Figure 20. Comparison of total vorticity (plan view just below the critical layer) at t = 45 
between LES and DNS with fine grid [10]. 
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Figure 21. Rope-like structure derived from wall-normal density gradient at f = 45 for LES 
with grid 24 x 12 x 96 compared with corresponding structure from DNS with fine grid [10]. 
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Figure 22. Comparison of spanwise component of vorticity at t = 50 between LES and DNS 
with fine grid [10]. 
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Figure 23. Comparison of spanwise component of vorticity at t = 53 between LES and DNS 
with fine grid [10]. 
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Figure 24. Comparison of spanwise component of vorticity at t = 54 between LES and DNS 
with fine grid [10], 
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Figure 25. Comparison of spanwise component of vorticity at t = 55 between LES and DNS 
with fine grid [10]. 




Figure 26. Structure of spanwise vorticity of flow field in transition region with LES and con- 
traction Cj. 
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